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The Rice-Sneddon model for BaBi03 is a nice model Hamiltonian for considering the properties 
of polarons and bipolarons in a three-dimensional oxide crystal. We use exact diagonalization 
methods on finite samples to study the stability and properties of polarons and bipolarons. Because 
polarons, when they form, turn out to be very well-localized, we are able to converge accurately 
our calculations for two-electron bipolaron wavefunctions, accounting for the Coulomb interaction 
without approximation. Some of our results are compared with and interpreted by reference to the 
variational method of Landau and Pekar. We calculate both electronic and vibrational excitations 
of the small polaron solutions, finding a single vibrational state localized with the full symmetry 
of the polaron, which has its energy significantly increased. Both on-site (Hubbard) and long- 
range Coulomb repulsion are included in the bipolaron calculation, but due to the high degree of 
localization, the long-range part has only a small influence. For a reasonable on-site repulsion U 
equal to 2 times the band width W , bipolaron formation is significantly suppressed; there is a large 
window of electron-phonon coupling where the polaron is stable but the bipolaron decays into two 
polarons. 



PACS numbers: 71.38.+i, 71.30.+h, 63.20.-e 
I. INTRODUCTION 

The electron localized on an ion (or a few ions) can 
cause displacements of neighboring ions from their po- 
sitions in the crystal lattice. The quasiparticle formed 
by an electron and corresponding lattice displacements 
is called a polaron (!]]. The polaron is a small polaron if 
the electron is localized on a few ions or a large polaron 
otherwise |2). It is possible that the surrounding ions 
"overscreen" the negative charge of the localized elec- 
tron and attract another electron to the same site. Two 
electrons bound in such way form a bipolaron The 
important issue for bipolaron formation is a quantitative 
analysis of the relative strength of the electron-lattice in- 
teraction responsible for two electrons being coupled and 
the electron-electron (Coulomb) forces which try to break 
the bipolaron apart. If the Coulomb repulsion between 
electrons exceeds some critical value the bipolaron is less 
stable than two separated polarons. 

In present paper we study the polaron and bipolaron 
formation in a model originally proposed by Rice and 
Sneddon j| for doped BaBiOs, such as Bai_ x K x Bi03 
(BKBO). These materials are nearly cubic perovskites, 
with a fairly simple set of valence electron states, ex- 
hibiting superconducting transition temperatures as high 
as 30K at optimal doping. We choose this model both 
because of its relation to physically important materials 
(SrTiC>3 and WO3 could be studied by a closely related 
model with T2 g (i-states instead of s-electron states) and 
also because it makes a nice test model for looking at 
the criteria for polaron formation and the properties of 
polarons. In this study, for simplicity, we consider the hy- 
pothetical case of an almost empty band, corresponding 
to KBi0 3 , or BKBO with x=l. 

Our principal results are (1) there is a sudden jump 



at a critical coupling strength from a large (delocalized) 
polaron to a very small, well-localized polaron; (2) the 
bipolaron breaks apart even at moderate strength of on- 
site Coulomb repulsion; and (3) small polaron formation 
is accompanied by a characteristic formation of localized 
vibrations of shifted energy which may serve as a good 
experimental signature. 



II. MODEL HAMILTONIAN 



The Rice-Sneddon Hamiltonian is 

H = H t + H e - p h + H p h + He 



(1) 



where H t , H e ^. p h, H p h and He are , correspondingly, 
hopping, electron-phonon, phonon, and Coulomb terms. 
For the hopping term, only one bismuth s-orbital per cell 
is used, with hopping only to nearest bismuth neighbors, 
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This has the usual cosine dispersion relation with band- 
width, W, equal to 12i. We take the value of t to be 
200 meV following band-structure calculations ||. We 
assume that the same model for electrons in BKBO can 
be applied to all the region of K concentration, because 
the antibonding Bi(6s)-0(2p) conduction band in cubic 
BKBO is minimally affected by substitutional K doping 
at the Ba sites ||. First-principles LMTO calculations 
for different phases of BKBO: KBi03, Ko.sBao.sBiOs and 
BaBi03 reveal a single band near Fermi energy to be 
largely independent on potassium doping ||. 

Hypothetical KBi03 in this model has no electrons. 
Each Bi atom is in the 5+ ionization state with no re- 
maining valence electrons, whereas BaBi03 in this model 
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has a half filled band of Bi 4+ ions. Real BaBiC>3 is insu- 
lating, and is most simply understood as having alternat- 
ing Bi 3+ and Bi 5+ ions. The Bi 3+ ions have two valence 
electrons bound to them, and repel the surrounding 6 
oxygen atoms, which are attracted to the more positive 
Bi 5+ ions. A nice term for this state is that it is a "bipo- 
laronic crystal." The Rice-Sneddon Hamiltonian incor- 
porates this possibility through the terms 
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The notation Ui yCC refers to the displacement of the oxy- 
gen atom (labeled by { i, a } ) which neighbors the «-th 
Bi atom in the a=x, y, z Cartesian direction. Only dis- 
placements in the direction a of the bond are considered, 
because these are expected to dominate the physics of 
polaron formation. The notation Ui + ^ a means exactly 
the same as Ui jQ , whereas Ui_ jQ refers to the displace- 
ment of the oxygen atom which neighbors the i-th Bi 
atom in the a = — x, — y, — z Cartesian direction. The 
atoms labeled by (i— , a) can also be labeled in the form 
(«', a) by reference to the appropriate nearby Bi atom 
i'. Eq. (||) contains the effect that an "inhaling" of the 
negative oxygen ions around a central Bi ion will raise 
the on-site energy of the bismuth s-orbital, by amount 
equal to g per fractional displacement u/a of each of the 
6 surrounding atoms. However, this costs elastic energy 
as given in Eq. (||). We take M as the oxygen atomic 
mass, and loq to have the value 65 meV of a typical oxy- 
gen bond-stretching vibration. (It corresponds to spring 
energy Mu^a 2 /2 being as 4 eV.) Various values of g will 
be used, in the physically expected range of 1-3 eV. The 
Bi-Bi interatomic distance a — 4. 28 A will be used as the 
unit of length, and the hopping parameter t is used as 
the unit of energy. Finally, there is a Coulomb interaction 
between electrons, 



^E 



E 

hi 



en 



(5) 



Various values of U comparable with W = 12t will be 
used for the on-site (Hubbard) term in the Coulomb in- 
teraction, and the value e = 5 will characterize the long 
range Coulomb repulsion. 

The vacuum of this model corresponds to KBi03, 
with no electrons, and only zero-point vibrational energy 
3Nhu>o/2. We will sec that shifts of the zero point energy 
are not very important, so we can ignore this term and 
define the vacuum as having energy zero. For the half- 
filled case, the ground state (for not too strong Coulomb 
repulsion compared to electron-phonon coupling) has a 



period-doubling distortion of the oxygens. Alternate bis- 
muth sites have the oxygens "breathing" either in or out, 
and electron charge either diminished or increased from 
the average of one electron per site. This corresponds to 
the experimental insulating state of BaBiOa, and gives 
a new vacuum into which carriers can be introduced by 
doping. Neglecting Coulomb repulsion, this regime was 
was studied numerically by Yu et al 0, who found, in 
agreement with experiment, that the insulating gap per- 
sisted for a wide range of K-doping. In d=2, this model 
was studied by Prelovsek et al. using an adiabatic treat- 
ment of the phonon degrees of freedom and the Hartree 
approximation for the Hubbard term |8J. 



III. POLARON 

Our only approximation (apart from finite size errors 
which are well-controlled) is the Born-Oppcnheimer (adi- 
abatic) treatment of the vibrations. Inserting one elec- 
tron into the empty-band vacuum, and letting the oxygen 
atoms have some fixed distortion pattern {u}, we look for 
the lowest energy one-electron state with wavefunction 
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where a,i({u}) are the site amplitudes of the electron wave 
function. Later the dependence of dj etc. on the param- 
eters {u} will be implicit and not explicitly designated. 
This electron state has energy eo({u}), measured relative 
to the bottom of the band, eo({0}). The total energy is 
this plus the elastic energy (H p h({u})). Then we vary 
the displacements {u} looking for the absolute minimum 
total energy. If the coupling constant g is small, the 
minimum occurs at {u} — {0} and has total energy 0. 
This corresponds to a large polaron solution, which in 
adiabatic approximation is just an electron in the bot- 
tom of the band of the undeformed crystal. If we were to 
include the non-adiabatic coupling of this electron with 
virtual phonons, there would be an alteration of the mass 
and energy of this electron. Specifically, the energy shift 
would be 
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in terms of the one electron energies ej and the phonon 
energies TitOk of the unperturbed band. This sum can be 
evaluated as follows, 



Ae 
t 



(f) ; 



Hujq 



2Muj%a 2 



l-S 



hujQ 

~4T 



N i—- 1 x 



f(Q)=sin 2 (- 



/(Q) 



^) + , sm 2 ( ^ ) + STO 2 ( 



Q z a 



(8) 
(9) 

) ■ (10) 



2 



4.0 



I 2.0 

(0 

DC 



0.0 




(a) 



>> 
cn 
i- 

c 

UJ 



24 



Q. 

(0 

C5 12 







(b) 


\ V 


(c) 


f 




5 10 15 

Cluster Radius 



20 



5 10 15 

Electron-Phonon Coupling 



FIG. 2. Tests of finite-size errors for different values of elec- 
tron-phonon coupling constants. The cluster radius is in the 
units of the Bi-Bi distance. The clusters used are No x No x JVo 
and (No — 1) x No x (No + 1), out to a maximum of iVrj of 20. 
The cluster radius is in the former case or N 1/3 (Nq - 1) 1/3 
in last case. 



FIG. 1. Calculated properties of a polaron: (a) radius 
rp/a, (b) total energy per electron E to t/t, and (c) the gap 
(ei — £2)/£ in the electron spectrum as a function of elec- 
tron-phonon coupling g/t. The circles give results for the 
single-electron polaron and the diamonds give the radius, half 
the total energy, and the electronic gap for the bipolaron. The 
Coulomb repulsion is omitted for the bipolaron (U = U = 0). 

For our choices of ui and t, the ratio Ucuo/At is 0.081 and 
the sum S in Eq. (Q) is 0.05. Thus the self-energy shift 
of the large polaron is w 4 x 10~ 3 x (g/t) 2 which turns 
out to be small compared to the energies that we will 
find for the small polaron regime. Thus we can safely 
ignore the non-adiabatic effects. By a similar argument 
(which wc will explain in more detail later) the zero-point 
contribution to the elastic energy can be ignored, and our 
elastic contribution to the small polaron energy is just the 
second term of Eq. (||). 

To evaluate the one-electron energy eo({u}) for the 
distorted lattice requires a finite size system, which we 
choose to be an orthorhombic cell (our "supercell" ) with 
N = Ni x iV2 x A*3 Bi atoms on a cubic lattice, and 3N 
oxygens on the Bi-Bi bonds, and periodic boundary con- 
ditions. The Lanczos technique |J was used for finding 
the ground state energy and a few lowest excited states of 



the Hamiltonian (|l|), and conjugate gradient minimiza- 
tion was used to find the optimum values of the oxygen 
displacements {u}. Beyond a critical value (<j>/i) c =8.57 
it becomes favorable for oxygens to distort and form a 
localized small polaron state. We define the location fo 
and radius rp of the polaron by 

i 

r|> = £k| 2 (^-n)) 2 . (12) 

i 

The radius of the polaron at the transition is 0.49 in 
units of Bi-Bi distance, that is, it is very well-localized 
with 90% of the electron density concentrated on one 
site. As g increases beyond the critical value, the radius 
further shrinks, and the binding energy rapidly increases 
to values of order t and bigger. 

Our results are plotted in Fig. [j]. For g/t less than 
the critical value, the radius is shown as a finite number, 
w 4a, reflecting the finite size of the cell; the actual radius 
is infinite. For values of g/t slightly less than critical, 
our minimization proceedure locates a metastable small 
polaron solution with a small positive energy, which is 
shown in Fig. [l] as a small hysteretic region. 
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FIG. 3. Total energy versus polaron radius calculated in 
Landau-Pekar approximation for electron-phonon coupling 
constants g/t—0., 8.23, 10., 12., and 15. from top to bottom. 
The energy minima for g/t= 10., 12., and 15. correspond to 
energies of stable localized states. The extremum at g/t=8.23 
corresponds to a metastable localized state. Results from fi- 
nite cluster diagonalization are shown as X's. For numerical 
comparison see Table 1. The inset expands the large rp part 
of the same curves. 



Because the small polaron is so well-localized, the error 
in our calculation due to the finite size supercell is easy 
to control. To test this, we have varied the size of the 
supercell from a minimum of 2 x 2 x 2 to maximum of 
20x21x22. The results are shown in Fig. g. The polaron 
radius and total energy are insensitive to cluster size if 
the number of Bi-atoms is > 200. Near the transition, 
for cluster size not too big, the transition onset varies 
with cluster size. The total energy always diminishes 
with increase of cluster size until it becomes independent 
of cluster size. At g/t far enough from the critical value 
the results are almost the same for all the clusters sizes. 
The results of Fig. [I] have no noticeable size dependence. 

The non-adiabatic corrections to the large polaron en- 
ergy are 0.0064i at the transition point, which gives an 
unimportant correction to the critical value of g/t. The 
small polaron solution is 2iV-fold degenerate: it can form 
at any of the N Bi sites, with either spin. In this paper we 
ignore another non-adiabatic effect, the weak vibration- 
assisted tunnelling which lifts the translational degener- 
acy to make a narrow band with only spin degeneracy 
remaining. 

We now compare our numerical results with analytic 
results obtained by a variational method introduced by 
Landau and Pekar (LP) ]lQ , |Tl|] . The electronic wave 
function is chosen to have Gaussian form 



TABLE I. Comparison of results obtained by Lan- 
dau-Pekar (LP) variational method and by exact diagonaliza- 
tion of finite clusters (cluster). The electron-phonon coupling 
constants are in units of bandwidth (W — 12t). The first 
row in the table corresponds to the critical coupling constant 
at which a metastable localized state occurs. The total ener- 
gies -Ej^tai an d ^totai" are m units of t. The polaron radius 
is given in units of the Bi-Bi distance. Cluster calculations 
give a slightly smaller total energy and larger radius. The 
phonon spectrum is perturbed by the transition from delocal- 
ized state to localized (see text). In LP approximation only 
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where Co is the normalization constant and (3 is the varia- 
tional parameter which we call the LP parameter. There 
are now two sequential minimizations to perform [p"2|| . 
First for fixed (3 the optimum displacements {u((3)} are 
found. Then these are used to evaluate the trial total 
energy E((3), and a second minimization is performed to 
find the optimum (3. We find analytic formulas for E((3) 
and the polaron radius rp((3), 
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(14) 
(15) 



where q = exp(— (3 2 ) and 9i{q) = 0^(0, q) are Jacobi's 
theta functions |L3|]. 9' 3 (q) is the derivative of 3 (q). 
These equations define an implicit function E(rp) which 
is plotted in Fig. ||for the values g/t equal to 0., 8.229, 
10., 12. and 15. X-marks indicate the exact results from 
our finite cluster calculations. The agreement is remark- 
able, especially, for the values of minimal total energy 
(see Table 1). The LP solution gives a smaller value 
of the polaron radius. The variational solution for an 
infinite system agrees with the exact solution on finite 
clusters in finding a first-order large to small polaron 
transition, with no intermediate regime of large polarons 
exist. Figure |^ explains the hysteresis found in the nu- 
merical results obtained by exact-diagonalization. For a 
small range of g just below the critical value the small 
polaron state is locally stable but separated by an energy 
barrier from the global (delocalized) minimum. The nu- 
merical solution follows this metastable branch until it 
disappears. 



ai = C (/3)e X p(-/3 2 (r l -f ) 2 /2) 



(13) 
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IV. ELECTRONIC AND VIBRATIONAL 
EXCITATIONS OF THE POLARON 



An advantage of the exact diagonalization method is 
that is enables an equally good and easy calculation of 
electronic excitations of the Franck-Condon type where 
the lattice distortion is frozen in place. We simply ex- 
amine the next higher-lying eigenstates without further 
alteration of the parameters {u}. In the range of pa- 
rameter space we have explored, we have not encoun- 
tered a second bound state in the polaronic well. The 
electronic spectrum has a gap, and the minimum energy 
electronic excitation is a delocalized state. The energy 
of this transition, denoted the "gap" energy, is plotted in 
Fig. |l}c. At the onset of polaron formation, the gap has 
a value 5.00t which increases rapidly for larger coupling 
constants. There is a small but noticeable finite size error 
in the gap calculation since the lowest electronic excited 
state is extended to infinity, but cut off at the supercell 
boundary in our work. 

When a small polaron is formed, the interaction be- 
tween the localized electron and the lattice vibrations can 
cause both a renormalization of the electron energy and 
of the phonon energy. Referring to Eq. (Q), it is clear that 
the gap, or minimum value of eo — e j m the denominator, 
makes a change in the electron self-energy shift relative 
to the one already calculated for the delocalized large po- 
laron, probably reducing the shift because of the larger 
denominator (although matrix element changes need to 
be considered also.) However, since the shift is certainly 
small compared to the gap itself (of order t), this effect 
can be neglected. 

A more interesting effect is the change in the local 
vibrations near the localized electron. If we know the 
one-electron energies Cj and the corresponding states |j) 
at the optimal set of displacements {uo}, then standard 
perturbation theory for small displacements around these 
displacements gives 
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This equation omits terms containing second derivatives 
of H because they are consistently omitted in our model. 

The LP approximation gives a particularly simple so- 
lution to this problem. Since we don't have a complete 
set of states in this approach, instead, we find the energy 
as a general function of the displacements {uq} and the 
LP parameter /3, 

E(u,(3) = \tf-A -u + D((3) -u + f(f3) , (17) 

where u is the 3iV-vector displacement, A is the bare 
force constant matrix (which is a constant times the unit 
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FIG. 4. Phonon frequency shifts, uj l vh /uio, at coupling con- 
stants g/t=8.23, 10., 12., and 15. for cluster size 5x6x7. 
The phonon modes (circles) are distributed from ojq to oj ma x, 
the last being well separated from others. Results from the 
Landau-Pekar approximation are shown by diamonds. 



matrix in our model), L is the force on the atoms caused 
by the localized electron, and / is the localized electron 
hopping energy. Expressions for L and / are easy to de- 
rive. Straightforward linear algebra leads to expressions 
for the optimum values {uq} and Pq. We then Taylor- 
expand Eq. ([Tt]) to second order for small deviations 
{Su} and Sp around the optimum values. Finally, for 
fixed deviations {Su} the optimum value of (5/3 is found. 
Inserting this into the Taylor expansion, the modified 
force constant matrix A' is found, 



A' = A 



L'L 



.t 



Lt • A- 1 ■ L" - f" 



(18) 



The primes on the right hand side of Eq. (]l8j) denote 
derivatives by (3. Note that the alteration of the force 
constant matrix in LP approximation is factorizable, and 
since A is proportional to the unit matrix, only one eigen- 
value is altered, the corresponding eigenvector being pro- 
portional to L' . The static displacements {uq} in LP 
approximation are given by — A~ l L((3q). Thus in LP ap- 
proximation, one vibrational eigenvector splits off from 
the degenerate frequency ujq, shifting to higher energy, 
and having an eigenvector proportional to the derivative 
of the static displacements by the LP parameter (3. The 
symmetry of this mode is identical to the symmetry of 
the static displacement. 

We have also made an exact calculation of the mod- 
ified vibrational spectrum using finite clusters, and the 
answers, shown in Fig. |], and also in Table 1, agree nicely 
with the LP approximation, but in addition to the one 
strongly altered frequency, a few other frequencies are 
pulled weakly above the unperturbed frequency ujq. 
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FIG. 5. Total energy (per electron) of the bipolaron for 
different values of Coulomb repulsion. The supercell size is 
5x6x7. At U = 2W for g/t between 8.6 and 11.3, the total 
energy corresponds to two small polarons separated as far as 
possible in the cell; because the cell is not infinitely large, 
there is Coulomb repulsion which raises the energy above the 
isolated polaron energy shown as the thin curve (identical to 
Fig. 0). 



Thus we expect that a characteristic signature of the 
small polaron state should be a localized vibrational 
mode whose symmetry copies that of the polaronic dis- 
tortion, that is, the symmetry is the same as the point 
symmetry in the crystal of the ion on which the polaron 
is centered (full cubic symmetry A\ in our case.) Such 
modes might be measureable by Raman scattering using 
a laser which is resonant with an electronic transition of 
the polaron. Also they might appear as side-bands on 
the electronic polaron absorption spectrum. 



V. BIPOLARON 

We now ask what happens in our model when a second 
electron is added. If we neglect the Coulomb interaction, 
the answer is that two spatially separated polarons are 
unstable relative to formation of a singlet bipolaron state 
in which both electrons are on the same site. If we al- 
low no further lattice relaxation beyond the single elec- 
tron polaron, then the energy of the bipolaron is already 
lower than two separated polarons because the (negative) 
electronic eigenvalue is doubled but the positive lattice 
strain energy is unchanged; additional lattice relaxation 
will occur only if it lowers the energy, and since there are 
now two electrons exerting each force g on the neighbor- 
ing oxygens, there will be additional relaxation. Results 
are shown in Fig. [I] where we plot the total energy per 



FIG. 6. Bipolaron radius as a function of g/t for different 
values of Coulomb repulsion. Vertical lines indicate stability 
limits where bipolarons (with U decay into isolated small 
polarons, or into isolated large polarons when U>W. At g/t 
less than the stability point, radii of metastable bipolaron 
states are shown. The domain of metastability is artificially 
enhanced by finite size effects. 



electron. The critical coupling for bipolaron formation 
is <7/t=6.08, significantly less than for polaron formation 
which starts at g/t=8.57. At onset of bipolaron forma- 
tion, the radius and the electronic gap (0.44a and 5.00i, 
respectively) are approximately the same as for polaron 
formation at its onset, but at equal values of g the bipo- 
laron is smaller and has a larger gap. 

Of course it is unrealistic to ignore the Coulomb repul- 
sion which will act in the direction of destabilizing the 
bipolaron. Our model permits us to make exact (finite 
size) calculations for the bipolaron by solving the appro- 
priate two-particle equation, that is, finding the exact 
two-particle wavefunction 



*({«}) = 2 au({«»c: 



(19) 



This calculation is of course far more demanding than the 
corresponding polaron case, Eq. (JsJ) , because it requires 
on each step of the minimization procedure finding the 
smallest eigenvalues of a N 2 x N 2 matrix rather than a 
N x N matrix. Our Lanzcos algorithm has allowed us 
to calculate for N as large as 5 x 6 x 7. Because the 
bipolaron turns out to be again well- localized, the finite 
system size does not cause a noticeable error. Also for 
the same reason, the long-range Coulomb repulsion is not 
very important. 

Results for the bipolaron radius, energy, and excitation 



gap are shown in Figs, 
strength U (see Eq. 



||, and 0. The onsite Coulomb 
fel)) is estimated as e 2 /2rBi, while 
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the long-range Coulomb interaction is characterized by 
a parameter we call U' — e 2 /ea. For rrji = 0.5 and 
e = 5 one finds U = 3.5 eV = 1.46W and W = 0.7 eV 
= 0.29IF, where the bandwidth W is 12t or w 2.4 eV. 
Our calculations are shown for the cases U/W= 1 and 2 
with E7'=0, and for U/W= 2 with U'= 0.2. Fig. | shows 
that on-site Coulomb repulsion with the small value of 
U = W destabilizes the bipolaron until g/t reaches 8.95, 
slightly beyond the stability point g/t — 8.57 for the 
polaron. A more physical choice of U = 2W prevents 
stable bipolaron formation until g/t « 12.8. Our numer- 
ical solutions show bipolarons winning over separated po- 
larons at somewhat smaller values of coupling constants 
(g/t = 8.60 for U = W and g/t = 11.3 for U = 2W .) 
This is a finite size error caused by the fact that in our 
cell, two separated polarons are only separated by half 
the cell diagonal, and have a small repulsion which artifi- 
cially raises their energy, favoring bipolarons for smaller 
coupling than would be the case in a very large cell. But 
since we have already an accurate answer for the energy 
of two well-separated polarons, our numerics tells us the 
true stability point. As is seen in Fig. M, there is also 
a small hysteretic regime indicating that the transition 
from separated polarons to bipolarons is first order. As 
is seen in Fig. ^[ adding the long-range Coulomb term 
(U' = 0.2U) causes a small additional postponement of 
bipolaron formation, but hardly alters the properties for 
values of g/t where it is stable. 

The radius shown in Fig. ^ is defined exactly as in 
Eq. ( fl2"| ) after defining the probability \a[\ 2 for finding 
an electron on site i in the usual way, 

ki 2 = £k,f- ( 2 °) 

3 

Fig. ^ shows that the radius of the bipolaron is increased 
by increasing the value of U, but decreased by increasing 
the value of U' . In all cases the bipolaron has a smaller 
radius than the polaron at the same value of g. 

An interesting feature shown in Fig. fj] is that the gap 
is larger at the onset of bipolaron formation in the pres- 
ence of the long-range Coulomb repulsion, presumably 
due to stronger localization of electrons. At fixed g/t the 
Coulomb forces reduce the gap with increasing U . 
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FIG. 7. The gap in the electron spectrum as a function 
of g/t. Due to finite size of cluster the gap is finite even in 
delocalized state. 

solution. This jump is not caused by finite-size errors and 
is present also in variational calculations using Landau- 
Pekar approximation. The total energy has hysteretic 
behavior with metastable states occuring near the criti- 
cal coupling constant. These metastable states could in 
principal be observed, for example, by tuning the cou- 
pling constant with applied pressure. A gap opens in 
the electron spectrum at the transition from delocalized 
to localized polaron states, and new localized vibrational 
states occur with energies increased above those of the 
undoped host. 
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VI. SUMMARY 

Bipolaron formation is strongly affected by Coulomb 
forces in a cubic perovskite lattice. Due to Coulomb re- 
pulsion between two electrons localized on the same site 
the onset of bipolaron formation can be postponed and 
polaron states are energetically favorable. The polarons 
and bipolarons formed in this lattice are small and exist 
only above some critical value of electron-phonon cou- 
pling. The transition from delocalized to localized po- 
laron state is discontinuous, with no intermediate- size 
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